RSQSim on TPV15 branched-fault benchmark
RMarkdown setup
library(knitr)
knitr::opts_chunk$set(cache = TRUE, cache.lazy = FALSE, warning=FALSE,
fig.width=6, fig.height=6, autodep=TRUE)Functions to setup runs
Function to build the fault geometry
buildFault = function(l=400, setoff=0, rundir="TPV15", run=1) {
# through fault
x0 = 0
y0 = 0
z0 = 0
L = 28e3
W = 15e3
w = l
nl = round(L/l)
nw = round(W/w)
strike = 0
dip = 90
rake = 0
ddot = 0
sectNum = 1
system2("/Users/dinger/RSQSim.github.master/BuildFaultModels/singlePlanarFault",
args=c(x0, y0, z0, L, W, nl, nw, strike, dip, rake, ddot, sectNum),
stdout=file.path(rundir, "fault1.flt"))
# branch fault
strike = 30
x0 = setoff * sin(strike*pi/180)
y0 = 16e3 + setoff * cos(strike*pi/180)
z0 = 0
L = 12e3
nl = round(L/l)
nw = round(W/w)
strike = 30
sectNum = 2
system2("/Users/dinger/RSQSim.github.master/BuildFaultModels/singlePlanarFault",
args=c(x0, y0, z0, L, W, nl, nw, strike, dip, rake, ddot, sectNum),
stdout=file.path(rundir, "fault2.flt"))
system2("cat", args=c(file.path(rundir, "fault1.flt"),
file.path(rundir, "fault2.flt")),
stdout=file.path(rundir, paste0("faults.flt")))
unlink(c(file.path(rundir, "fault1.flt"),
file.path(rundir, "fault2.flt")))
flt = readFault(file.path(rundir, "faults.flt"))
return(flt)
}Function to set up initial stresses
Note that the 3 km x 3 km nucleation region centered on \(y\) = 8 km, \(z\) = -7.5 km may or may not fit exactly in a mesh with a given side length.
initStresses = function(flt, sigma=rep(120, length(unique(flt$iseg))),
tau=rep(70, length(unique(flt$iseg))),
ynuc=8000, znuc=-7500, wnuc=1500, segnuc=1, taunuc=81.6) {
sigma.out = sigma[flt$iseg]
write(sigma.out, file=file.path(rundir, "initSigma.txt"), ncol=1)
tau.out = tau[flt$iseg]
nuc = which(flt$iseg == segnuc & abs(flt$z - znuc) <= wnuc & abs(flt$y - ynuc) <= wnuc)
tau.out[nuc] = taunuc
write(tau.out, file=file.path(rundir, "initTau.txt"), ncol=1)
}Elastic modulii
These are simple and the same for all the runs we’ll be doing here.
TPV15 (and TPV15) use \(\rho = 2670 \mathrm{kg}/\mathrm{m}^3\), \(V_s\) = 3464 m/s, and \(V_p\) = 6000 m/s. RSQSim needs the Lame parameters \(\lambda\) and \(\mu\). In terms of the density and wave speeds these are
\[\begin{align} \mu &= \rho V_s^2 \\ \lambda &= \rho \left( V_p^2 - 2V_s^2 \right) \end{align}\]rho = 2670
Vs = 3464
Vp = 6000
lameMu = rho * Vs*Vs / 1e6 # in MPa
lameLambda = rho * (Vp*Vp - 2*Vs*Vs) / 1e6Friction parameters and initial state
We will use our usual fixed values of \(V^\mathrm{EQ}\) (though will definitely want to try variable \(V^\mathrm{EQ}\) eventually) of 1 m/s (update: given the specified \(V_s\), \(\Delta\tau\), and shear modulus, the shear impedance relation implies a slip speed of 1.57 m/s, but I have adjusted it to 1.2 m/s to better match the rupture propagation speed seen in the dynamic models), \(V^\ast\) of \(10^{-6}\) m/s, \(D_c\) of \(10^{-5}\) m, \(a\) of 0.01, and \(b\) of 0.015. Then we will choose \(\mu_0\), and \(\theta_0\) from the two equations in http://rsqsim.ucr.edu/wiki/index.php/Stress_drops_in_artificially_nucleated_events#Further_attempt_at_more_accurate_estimate_of:
\[\begin{align} \mu_0 &= \frac{\tau_f}{\sigma} + \left(b-a\right)\ln\left(\frac{V^\mathrm{EQ}}{V^\ast}\right) \\ \theta_0 &= \left(\frac{D_c}{V^\ast}\right)^\frac{b-a}{b-a+a^\prime} \left(\frac{l b \sigma}{G V^\mathrm{EQ}}\right)^\frac{a^\prime}{b-a+a^\prime} \exp\left(\frac{\tau_p/\sigma - \mu_0}{b -a +a^\prime}\right) \end{align}\]And then force the nucleation region to fail immediately at its initial shear stress level by lowering \(\theta\) there by the appropriate amount. RSQSim can do that automatically, but I’ll just calculate the numbers here and write out an initial \(\theta\) file instead.
The approx for \(\tau_f\) is known to be pretty rough, so after an initial run showed that the final \(\tau\) was about 2.6 MPa too high, I adjusted the \(\tau_f\) in the formula for \(\mu_0\) above by this amount (in code block below).
For TPV15 (and TPV15) we have that \(\tau_p = \mu_s \sigma = 0.677 \cdot 120 \mbox{ MPa} = 81.24 \mbox{ MPa}\) and \(\tau_f = \mu_d \sigma = 0.525 \cdot 120 \mbox{ MPa} = 63 \mbox{ MPa}\).
fricTheta0 = function(flt, rundir, mus=0.677, mud=0.525, sigma=120,
VEQ=1.22, Vstar=1e-6, Dc=1e-5, a=0.01, fa=0.1, b=0.015,
dtauf = -2.6,
ynuc=8000, znuc=-7500, wnuc=1500, segnuc=1) {
taup = mus * sigma[1]
tauf = mud * sigma[1]
aprime = fa*a
mu0 = (tauf + dtauf)/sigma[1] + (b-a)*log(VEQ/Vstar)
theta0 = (Dc/Vstar)^((b-a)/(b-a+aprime)) * ((l*b*sigma[1])/(lameMu*VEQ))^((aprime)/(b-a+aprime)) *
exp( (taup/sigma[1] - mu0)/(b - a + aprime) )
theta.nuc = RStheta(Vstar, mu0, taup, sigma[1], a, b, VEQ, Dc)
theta = rep(theta0, flt$np)
nuc = which(flt$iseg == segnuc & abs(flt$z - znuc) <= wnuc & abs(flt$y - ynuc) <= wnuc)
theta[nuc] = theta.nuc
write(theta, file=file.path(rundir, "initTheta.txt"))
return(mu0)
}Building the RSQSim input file
buildRSQSimInputFile =
function(rundir, run, a, b, Dc, mu0, ddotStar, VEQ, lameLambda, lameMu) {
unlink("eqs..out")
system2("/Users/dinger/RSQSim.github.master/runRSQSim", stdout=NULL, stderr=NULL)
eqs = readEqs("eqs..out", saveRData=FALSE)
unlink(c("eqs..out", "eqs..out.RData", ".Screen.Time.out", ".times_saved.out"))
params = eqs$params
rm(eqs)
params$A_1 = a
params$B_1 = b
params$Dc_1 = Dc
params$mu0_1 = mu0
params$ddotStar_1 = Vstar
params$alpha_1 = 0.25
params$ddotEQ_1 = VEQ
params$lameLambda = lameLambda
params$lameMu = lameMu
params$nEq = 1
params$faultFname = "faults.flt"
params$outFnameInfix = paste0("TPV15.run", run)
params$writeTau = 2
params$writeSigma = 2
params$initTauFname = "initTau.txt"
params$initSigmaFname = "initSigma.txt"
params$initThetaFname = "initTheta.txt"
writeRSQSimInputFile(params, file.path(rundir, paste0("TPV15.run", run, ".in")))
}Runs with 400 m elements, VEQ = 1.22 m/s, and varying offsets of branch fault from main fault
Note that with 1 km elements and set up as below, rupture would not propagate away from nucleation region. Currently using 400 m elements.
Zero offset
run = 1
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
l = 400
setoff = 0
ynuc = 8000
znuc = -7500
wnuc = 1500
segnuc = 1
taunuc = 81.6
mus = 0.677
mud = 0.525
initSigma = c(120, 120)
initTau = c(70, 78)
VEQ = 1.22 # note: usual default is 1.0, shear impedance in this case says 1.57
Vstar = 1e-6
Dc = 1e-5
a = 0.01
fa = 0.1
b = 0.015
dtauf = -2.6
eqs = NULL
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)hist(eqs[[1]]$taupList, breaks="FD")
c(mean(eqs[[1]]$taupList), mus*initSigma[1])## [1] 81.47452 81.24000tauf.rsq = readFaultVal(eqs[[1]], val="tau", patches=sort(eqs[[1]]$pList),
writeVal="final")$val[1,]
hist(tauf.rsq, breaks="FD")
c(mean(tauf.rsq), mud*initSigma[1])## [1] 63.40268 63.00000100 m offset
run = 2
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 100
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)200 m offset
run = 3
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 200
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)300 m offset
run = 4
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 300
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)400 m offset
Only continued to this offset because this is one of our element sizes, which what TPV15 did - but with smaller elements
run = 5
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 400
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)Runs with 200 m elements, VEQ = 1.22 m/s, and varying offsets of branch fault from main fault
Zero offset
run = 6
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
l = 200
setoff = 0
VEQ = 1.22 # note: usual default is 1.0, shear impedance in this case says 1.57
b = 0.015
dtauf = -2.6
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)hist(eqs[[run]]$taupList, breaks="FD")
c(mean(eqs[[run]]$taupList), mus*initSigma[1])## [1] 81.56892 81.24000tauf.rsq = readFaultVal(eqs[[run]], val="tau", patches=sort(eqs[[run]]$pList),
writeVal="final")$val[1,]
hist(tauf.rsq, breaks="FD")
c(mean(tauf.rsq), mud*initSigma[1])## [1] 63.7541 63.0000100 m offset
run = 7
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 100
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)200 m offset
run = 8
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 200
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)300 m offset
run = 9
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 300
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)400 m offset
run = 10
rundir = file.path("TPV15", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 400
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV15.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV15.run", run, ".out"))
setwd(here)Rupture time plots
400 m elements
Zero offset
100 m offset
200 m offset
300 m offset
400 m offset
200 m elements
Zero offset
100 m offset
200 m offset
300 m offset
400 m offset
To do
- differing resolutions
- variable slip speed